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, In this paper we present an approach aimed at performing many-body calculations of Born- 

effective charges of crystalline insulators, by including the electron-correlation effects. The scheme 
is implemented entirely in the real space, using Wannier-functions as single-particle orbitals. Corre- 
lation effects are computed by including virtual excitations from the Hartree-Fock mean field, and 
the excitations are organized as per a Bethe-Goldstone-like many-body hierarchy. The results of 
our calculations suggest that the approach presented here is promising. 
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The Born effective charge (BEC) of a periodic solid is an important phenomenological quantity which connects the 
electronic structure of the system to its phononic properties*!- Of late, in the context of ferroelectric materials and their 
phase transitions, BEC has generated tremendous amount of interest^ Using BEC, one can also describe the lattice 
dynamics, and its coupling to infrared radiation, in a simple intuitive manner^ Most of the modern calculations 
of BECs are based upon the Berry-phase-based theory of macroscopic polarization developed by King-Smith and 
VanderbiltjA The aforesaid formalism is based upon single-particle orbitals, and therefore, can be implemented in a 
straight-forward manner within ab initio density-functional theory (DFT) fi£ or the Hartree-Fock (HF) framework^ As 
far as many-body calculations of polarization properties are concerned, Martin and coworkers have proposed several 
approaches which, to the best of our knowledge, have not been implemented within an ab initio methodologyAi 
Filippetti and Spaldin have recently implemented an ab initio method aimed at including correlation effects by using 
'"O ■ a self-interaction-corrected (SIC) density-functional approach^. 

Recently, we have developed a wave-function-based ab initio methodology aimed at performing electronic structure 
calculations on crystalline insulators i^dft' 11 ' 12 The approach uses Wannier-functions as single particle orbitals obtained 
at the Hartree-Fock level, which can subsequently be used to include the electron correlation effects, if needed. The 
approach has been applied to calculate ground state geometries, cohesive energies, and elastic properties of crystalline 
insulators at the Hartree-Fock leveljSiifi as well as at the correlated levelpiLA^ Moreover, within the Berry-phase 
formalism of King-Smith and Vanderbilt^ we have also used our approach to compute the BECs of several ionic 
insulators at the HF levels The purpose behind the present work is to use our Wannier-function-based methodology 
' to perform correlated calculations of the BECs of insulators. Since, ours is a real-space approach, we start with the 
following expression for the electronic contribution to the polarization per unit cell (P e ) valid for insulators^ 
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where A is a parameter governing the state of crystal (for the present case, it represents atomic displacements), Q 
is the volume of the unit cell, q e is the electronic charge, N (— > oo), represents the total number of unit cells in the 
crystal, R e = J2k=i rk ^ s tne many-particle position operator for the N e electrons of the crystal, and |\&q^) represents 
the correlated ground-state wave-function of the infinite solid. Next we verify that for an infinite crystal, Eq. (TIJ 
above is consistent with the Berry-phase-based expression for the BECs derived by King-Smith and Vanderbilfe^ 
at the mean-field level. If we express the mean-field (HF or otherwise) ground-state many-particle wave function 
|$q )of a crystal in terms of Wannier functions expressed in terms of square- integrable occupied Wannier functions 
{W„(r — Rj), i = 1, . . . , N, n = 1, . . . , M} located in the N unit cells constituting the solid^fi then using the Slater- 
Condon rules governing the matrix elements of a one-body operator between two many-particle states^ we obtain 
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($i x) \R e \4 x) )=Nj2fn f v\W n (r)\ 2 dr + n e (f2^i), (2) 

n=l J i=l 

where {Ri, i = 1, . . . , N} are the lattice vectors of N unit cells of the crystal, W n (r) is the n-th Wannier function of 
the unit cell, /„ is the number of electrons in the n-th Wannier function (/„ = 2, for band insulators), M is the total 
number of occupied Wannier functions per unit cell, and n e — 53 n =i ls ^ ne total number of electrons per unit cell. 
However, if the N unit cells are distributed among complete shells (stars), then 53i=i ^ = 0- Combining this result 
with Eqs. (01 and iTjQl, we obtain the mean-field expression (P^) for the polarization per cell for a crystal 
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P {X) = qe/n^fn / r|W^(r)| 2 dr. (3) 

n=l J 

This equation is nothing but the Wannier-function version of the Berry-phase-based (mean-field) expression for macro- 
scopic polarization derived by King-Smith and VanderbiltjA who gave it an intuitive interpretation as being a sum 
over centers of Wannier functions of the unit cell. Note that expressions above are valid only in a real-space based 
approach where square-integrable Wannier functions are used as single-particle orbitals. If one were to use Bloch or- 
bitals instead, the expectation value of position operator will have to be computed differently^ Having demonstrated 
the equivalence of our starting expression (Eq. QJ) to the traditional theories at the mean-field level, we next examine 
its implications when a many-body expression for the ground state wave function IvE^), expressed in terms of virtual 
excitations from the mean-field wave-function, is usecU^ 

= C&W>+ E C$&|«W(n->a> + 

E C { ^ a ^\n^a- m ^(3)) + --- , (4) 

where the Greek indices a,/3,... represent the virtual Wannier functions while the Latin indices m, n, . . . represent 
the occupied ones. \&Q X \n —> a)) denotes a singly-excited configuration obtained by promoting one electron from 

the occupied Wannier function labeled n, to the virtual Wannier function labeled a. Similarly, |$q (n — * a \ 171 ~~ * 
(3)) represents a doubly-excited configuration with electrons being promoted from Wannier functions m, n to a, (3. 
Noteworthy point is that the occupied (m, n, . . .) and the virtual Wannier functions (a, 0,...) could be located in any 
of the unit cells of the infinite solid. The coefficients {Cq X \ C^a, C^ n . a p, . . .} can, in principle, be obtained using 
various available many-body techniques such as the configurations-interaction (CI) method, perturbation theory, etc. 
Next, we examine, the nature of contributions to the polarization vector P^ A ) arising from virtual excitations when 
a correlated wave function ( | q A ^ ) ) of the type of Eq. Q is used in Eq. QJ. In order to simplify things, we restrict 

our discussion to the contribution of the singly-excited configurations |$q \n — ► a)), although in our calculations 
all possible excitations needed to compute both one- and two-body increments have been considered (see discussion 
below). Thus, the expectation value of the dipole operator for a singly excited many-body wave function (assuming 

that |*[, A) ) is real) is 

(^ A) |R e |v|/( A )) = A^ A)2 (r)< A) + 

2V2NY / C£lC { X) (a\r\n(o))W +■■■, (5) 

where (r)^ = X^fiLi I r \W n (r)\ 2 dr (the HF expectation value per unit cell), and, using the translational symmetry, 
the sum over occupied Wannier functions has been restricted to those in reference unit cell, denoted as \n(o)), while the 
virtual Wannier function \a) can be in any unit cell of the solid. It is clear from Eq. (JJ that: (a) the expectation value 

of the dipole operator of the entire solid (\E / o A ^|R e | , I'o A ^) scales linearly with N as it should, and (b) the correlation 
corrections to the dipole moment/cell such as the second term of Eq. can be seen as due to the interactions 
between the electrons of the reference unit cell with those in the rest of the solid. Although, Eq. iJ^J has been derived 
for correlated wave functions containing only singly-excited configurations, however, it is easy to verify that even for 
more complex wave functions, only two other types of matrix elements, viz., (a|r|/3(o)) and (m\r\n(o)) contribute to 
the correlated expectation value. These matrix elements originate from interaction among different types of excited 
configurations. Because of the localized nature of the orbitals used, the dipole matrix elements will fall to zero rapidly 
with the increasing distance between the orbitals involved. For example, for a typical insulating solid, the dipole 
matrix elements are negligible when the orbitals involved are farther than third-nearest neighbors. 

Calculation of the correlated many-body wave function (Eq.QJ) °f an infinite solid is an extremely difficult task, 
thereby rendering the direct use of Eq. Q even more cumbersome. However, Stolid proposed the use of an "in- 
cremental" method of calculating correlated total energy (and wave function) of extended systems based upon a 
Bethe-Goldstone like expansion of correlation contributions. The approach was subsequently implemented by us to 
the case of infinite systems, and utilized to compute the total energy/cell and related quantities of bulk LiH^ii and 
several polymers^ In the approach, the total energy/cell is written as E ce u — Ehf + E corr , where Ehf is the HF 
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energy/cell of the system, and E corr is the contribution of correlation effects to the total energy/cell, computed as 

Ecorr ^ ^ ~t~ ^ ^ ^^ij 

i ' i^j 

| E Ae yfe + --- (6) 

where e^, Aey, Ae^fe, . . . etc. are respectively the one-, two- and three-body, . . . correlation increments obtained by 
considering simultaneous virtual excitations from one, two, or three occupied Wannier functions, and i, j, k, . . . label 
the Wannier functions involved^ However, using the incrementally computed many-particle wave function to compute 
the expectation value in Eq. Q is a tedious task which we avoid by using generalized Hellman-Feynman theorem, and 
the finite-field approach to compute dipole expectation values^ Accordingly, we perform the incremental calculations 
of energy/cell with the modified Hamiltonian 

H'(£,X) = H (X)-q e £-J2r k , (7) 

fc=i 

where -Ho (A) is the usual Born-Oppenheimer Hamiltonian for the solid with the given value of A, and £ is a user 
specified external electric fields From E^L(£) so computed, one can easily obtain 

N - d£~> (8) 

where j represents the Cartesian spatial component (j — 1,2,3). The partial derivatives in Eq. (jHJ, were computed 

numerically by performing the calculation of E^fa (£) for several small values of electric field £. The Wannier functions 
used in the present work were obtained by solving HF equations in the presence of an electric field, and thus are different 
from the ones used in our earlier works. 

Next we present and discuss our results for the cases of bulk LiH and LiF. In the present study also we have used 
the lobe-type contracted Gaussian basis functions used in our earlier worksi Si 1 ^ 11 ! 12 ! 1 ^ Unit cell Wannier functions 
for both the materials were described using basis functions centered in cells as far as the third-nearest neighbors of 
the reference celliSiifl For LiH we performed the calculations using the optimized lattice constant 4.067 A obtained 
in our earlier correlated calculation^ which is in excellent agreement with the experimental value of 4.06 A. For LiF 
we used the experimental lattice constant of 3.99 A. For both the systems, fee geometry was assumed, and anion 
and cation locations in the primitive cell were taken to be (0,0,0) and (a/2,0,0), respectively, where a is the lattice 
constant. For LiH, in correlated calculations Is Wannier function localized on Li as was treated as core and was 
held frozen, while for LiF, Wannier functions corresponding to Is orbitals of both Li and F were frozen during the 
correlated calculations. For computing BEC's, the parameter A corresponds to atomic displacements Au which was 
taken to be 0.01a (a is the lattice constant) in the x direction for the anion. The BEC tensor for cubic materials is 
diagonal and for the z-th atom of the cell it has only one unique component say Z*{i). It was computed using the 
formula Z*(i) = Zi + (Sl/e)-^, where Zi is the nuclear charge of the displaced atom, e = \q e \, and AP is the change 
in the unit cell polarization due to the atomic displacement Au. For Hellman-Feynman calculations of the dipole 
expectation value, we used the central difference formula, with the values for the external electric field £ — ±0.001 
a.u. in the x direction. 

First we verify whether the BEC's computed using the Wannier-function centers (Eq. Q) agree with those computed 
using the Hellman-Feynman theorem (Eq. <|HJl), at the HF level. Good agreement between the two calculations at 
the HF level will be a vindication of our approach, while any serious disagreement between the two results will be 
a setback, and will render further correlated calculations meaningless. The results from the two calculations are 
presented in tableland it is clear that the values obtained by the two methods are in excellent agreement with each 
other. Next we present the results of our correlated calculations in table |H| The table presents the changes in the 
values of BEC's as correlation effects of increasing complexity are included using the aforesaid incremental scheme. 
The many-body approach used to compute various correlation increments of Eq. (O was the full-CI method, as in 
our earlier worksi 11 ^ 2 In case of LiH we have performed correlated calculation including up to third-nearest neighbor 
(3NN) two-body correlation effects (ii while for the case of LiF these calculations have been restricted to the nearest 
neighbors (NN) two-body increments. The reason behind restricting the correlation effects for LiF to NN pairs is 
because the contributions beyond that (2NN, 3NN, ...) were found to be negligible. This is due to the fact that for 
LiH, the hydrogen anion is more diffuse as compared to the fluorine anion of LiF. Thus the valence electrons of LiH 
are comparatively more delocalized as compared to those of LiF, thereby making the correlation effects relatively 
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Table I: Comparison of Hartree-Fock Born charges of Li computed using the Wannier center approach (cf. Eq. ), and via 
the use of Hellman-Feynman Theorem (Eq. 0). 



System 


Born Charge 


Wannier-Center Approach 


Hellman-Feynman Approach 


LiH 


1.0417 


1.0418 


LiF 


0.9986 


0.9983 



Table II: Influence of electron correlation effects on the Born charge. Column with heading HF refers to results obtained at the 
Hartree-Fock level. Heading "one-body" refers to results obtained after including the corrections due to "one-body" excitations 
from each Wannier function of the unit cell, to the HF value. Two-body (O) implies results include additional corrections due 
to simultaneous excitations from two distinct Wannier functions located on the anion in the reference unit cell. Two-body 
(NN), two- body (2NN), and two-body (3NN) correspond to two-body correlation effects involving 1st, 2nd, and 3rd- nearest 
neighboring Wannier functions, respectively. 



System 


Born Charge 


HF 


one-body 


two-body (O) 


two-body (NN) 


two-body (2NN) 


two-body (3NN) 


Exp. 


LiH 


1.0418 


1.0302 




1.0193 


1.0183 


1.0003 


0.991" 


LiF 


0.9983 


0.9913 


0.9847 


1.0237 






1.045 6 



"Brodsky and Burstein— 
'Computed from the experimental data reported in Ref42. 



longer range in LiH. Inspection of table |H] reveals that for the case of LiH, at the HF level the BEC is overestimated, 
while for LiF it is underestimated. When the correlation effects within the reference unit cell are included for either 
of the systems, the value of the BEC decreases as compared to its HF value. This reduction can be seen as due to the 
mixing of the occupied Wannier function of anion with the unoccupied ones of the nearest-neighbor cation, termed 
"ion-softening", by Harrison^* However, the noteworthy point is that the ion-softening in the present case is being 
driven by the electron-correlation effects. As far as the longer range correlation effects (INN, 2NN, ...) are concerned, 
no clear trends are visible in table |n] For LiH we see monotonic decrease in the value of BEC with longer range 
correlation effects, while for the case of LiF the BEC increases as the nearest-neighbor correlation effects are included. 
However, in both the cases, upon truncation of the correlation series, the values of BEC's obtained are in excellent 
agreement with the experimental values. Finally, in order get a feel for the magnitude of correlation effects with and 
without the electric field, we present the values of various contributions to the correlation energy/cell for the two 
systems in table IIIII computed for the undistorted unit cell (A = 0). Results of similar calculations performed for the 
distorted unit cell are not presented here for the sake of brevity. From the table it is obvious that, as expected, the 
most important corrections to E ce u, due to the electric field are at the HF level. The correlation energies in nonzero 
field are reduced by small amounts as compared to their zero-field counterparts. Of course, these small changes in 
correlation energies for the A = and A ^ case in the end lead to the correlation corrections to the BEC's as 
depicted in table ITT1 

In conclusion, we have presented an ab initio Wannier-function-based many-body approach aimed at computing 
the Born effective charges of insulators. However, it is clear from the approach that it can also be used to compute 
other properties such as high-frequency dielectric constant, piezoelectric tensor, etc. of insulators. Work along these 
lines is presently under way in our group, and the results will be presented in future publications. 
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Department of Science and Technology (DST), Government of India. 
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£ 


Energy per cell (a.u.) 




(a.u.) 


Ehf 
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